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Abstract 

The aim of this article is to show, how computer algebra can be used when applying Liu's 
procedure. 

Although Mathematica® (a commercial product by Wolfram Research Inc.) is used in 
this article, it is possible to use other computer algebra systems as well. 

1 Thermodynamics 

1.1 Liu's Procedure 

(A very short introduction) 

Here a very short introduction will be presented, so that the reader will get a glance of what 
the Liu procedure is. For a detailed discussion of Liu's procedure and balance equations please 
consult the following literature: 0, [2j, 0, 0,0, 0, 



Balance equations 




Constitutive 
equations 








State Space 



Liu -Procedure 



Restrictions by 2nd Law 



Coleman-Mizel- 
Formulation of the 
2nd Law 



AH solutions of the 
Balance equations 
satisfy the 2nd Law 



Amendment to 
the 2nd Law 



There are no reversible 
process-directions in 
non-equilibrium 



Proposition 1 (Coleman-Mizel formulation of the 2nd Law) 

If Z is no trap, then the following statement is valid for all X : 



A(Z) • X = -B(Z) 
^a(Z)-X > -/3(Z) 



(1) 
(2) 



For all X, which satisfy the balance equations, the dissipation inequality is fullfilled automati- 
cally. 

Proposition 2 (Liu proposition) 0]/ 

Starting with the Coleman-Mizel formulation of the Second Law one can show that in large state 
spaces there exist state functions so that the following relations are valid. 

A(Z)-A{Z) = a{Z) 
-A(Z)-B(Z) > -P(Z) 

This means that there exist restrictions to the constitutive equations. 

1.2 Example for the Application of Liu's Procedure 
(for arbitrary non-relativist ic balances) 

Here a short example for general balances and a state space is given jH]. This example shall 
demonstrate how the balances on the state space and the matrix formulation is achieved. 
Therefore two balances and an inequality for arbitrary quantities are considered. 
The balances read: 

d t ai + V • (vai + b t ) = ci (3) 
d t a 2 + V • (va 2 + b 2 ) = c 2 (4) 
atf + V-(vf + 7) > (5) 

The state space is 

Z = {a u a 2 ,v,d t a 1 } (6) 
This is the domain of for the constitutive functions bi,b 2 ,^. 

Balances on the state space are obtained by using the chain rule when performing the deriva- 
tives. 

d t a 1 + V-(va 1 + b 1 (Z))-c 1 = (7) 

d t ai + V ■ Mi + v ■ Vai + 
dbi „ dbi „ dbi „ dbi „„ 

oai oa 2 ov oo t ai 

d t a 2 + V-(va 2 + b 2 {Z))-c 2 = (9) 
d t a 2 + V • va 2 + v ■ Va 2 + 

+ *. Vai + fi.Va 2 + fi.V 1!+ i^.Va iai -c 2 = (10) 

oai oa 2 ov oot a i 

d£(Z)+V-(v£(Z)+-y(Z)) > (11) 



dai da 2 dv dd t a\ 

+ V • v£ + u ■ t— Vai + u • Va 2 + «- 7r-% + «- ttt; V<9 t ai + 

aai oa 2 ov ao t a\ 

+^ L -Va 1 + p^-Va 2 + ^-Vv + ^-Vd t a 1 > (12) 

oai oa 2 ov oo t ai 
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As d t a>i is included in the state space it is not a higher derivative. 



Next the balances can be written in matrix formulation. 

A{Z)-X + B(Z) = 
a(Z)-X + (3(Z) > 

((X)) T = ( d t a 2 d t v d t d t ai Va x Va 2 Vv Vd t ai ) 
/ v + ^ || a4+^ 

((>*)) = 

\1 || , + H a 2 l + § $l.J 
((5)) = (^ Cl ) 

( ( n/)\ = ( M. <£ 9 Z v . M- + JtL v .J!L + JtL n +v .9i + <tL 

\ d(i2 dv ddtai — dai dai — da,2 da-2 s = — dv dv 

(09)) = (Jl-^x) 

Next the Liu equations can be written down 

MZ)-£{Z) = a(Z) 
-A(Z)-B(Z) > -(3{Z) 



A 2 





da 2 
dv 

dd t cii 



dbA fdb 2 \ d 1 

Ai v + - + A 2 - = v ■ — ^ + 



<9fli y \0<i\J dai da\ 

. dh . ( db 2 \ d£ d 1 

Ai- +A 3 V + 7 1 = + 



<9a 2 \ da 2 J da 2 da 2 

db A k ( db 2 \ ^ di d 1 

Ar^+A^ = v_-J^ + 



dd t ai dd t ai dd t a x dd t a\ 

The residual inequality reads 

di 

A x (a t ai - ci) - A 2 c 2 < ^-"^ a i 



Because of the principle of material frame indifference all derivatives with respect to v shall be 
set equal to zero: ^ = 0. The material frame indifference means that material laws shall be 
independent of the velocity of the observer. Using this the Liu equations read: 

° = I < 31 > 

o = T^r- ( 32 ) 

od t ai 

Al ^ + ^J +A2 ^J = ^da-^da-, (33) 

A 9lh A ( dl>2 \ v ^ (341 
1 da 2 2 \ da 2 J 

Ai(aii)+A 2 (a 2 i) = £ (35) 

<%i , <% 2 d£ i \ 

Ai— ^ + A 2 — ^ = <;._^ + _J_ 36 

The residual inequality remains unchanged 

Ai(ftoi - ci) - A 2 c 2 < — -9 t ai (37) 

aai 

From these equations the Lagrange parameters A(i£) can be determined. 

A 2 = |^ (38) 

aa 2 

da 2 

After inserting the Lagrange parameters into Liu's equations the restrictions to constitutive 
equations are given 



= I < 40 > 

= sab (41) 



e _^( ao l)V«il) _1 — + — ^ + — ^ = v^- + ^- (43) 

<9a 2 ^ =' J ^ =' da 2 da 2 \ da 2 J 

d£ ( A ( v-i ^ 3£ a& 2 ^7 

-^A))( a 'A) a^ + a^7 = fi \5^ + «Sr (44) 







The entropy production is given by the residual inequality 



This example was meant to demonstrate how the Liu equations and the residual inequality are 
achieved. This example shows also that the entropy production is independent of the highest 
derivatives included in the state space (this is a general property of the entropy density; see 
also [ID]). 



2 Application of Liu's Procedure using Mathematica 

For the correct handling of indices I will use the package "Ricci.m" by John M. Lee, which can 
be found on the MathSource server. 

«Ricci' 

2.1 Declaration of Tensors and State Space 

First it is necessary to define a bundle and declare the letters used for indices. 
Def ineBundle [f iber2,4,g, {i , j ,k, 1 ,m,n, o ,r ,y}] 

Next the tensors have to be defined. First I will declare the state space variables (zO, . . . , z5). 

Z = {nO,u,p,e, spa, spv} 



(* zO=nO *) Def ineTensor [n0,0] (* particle number density *) 
(* zl =u *) Def ineTensor [u, 1] (* four velocity *) 
(* z2=p *) Def ineTensor [p, 1] (* momentum flux *) 
(* z3=e *) Def ineTensor [e , 0] (* energy *) 

(* z4=spa *) Def ineTensor [spa, 2, Symmetries -> Skew] (* spin density *) 
(* z5=spv *) Def ineTensor [spv, 1] (* vector spin density *) 

Then the constitutive functions are declared: 

Def ineTensor [t , 2] (* stress tensor *) 

Def ineTensor [q, 1] (* heat flux *) 

Def ineTensor [scs, 3] (* couple stress *) 

Def ineTensor [spb, 2] (* 2-couple stress *) 

Def ineTensor [s , 0] (* entropy *) 

Def ineTensor [se , 1] (* entropy flux *) 

Def ineTensor [f , 1] (* external force *) 
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2.2 Balance Equations 



The next step is to write down the balance equations, which are for technical reasons resolved 
to zero, and the right hand side (= 0) is not written down: 
Balance of particle number density: 

d l (nu l ) = 



D[nO[x] u[U[l]] [x] ,x]) 
Balance of energy momentum 

d l {? l + \p i u l + \u i q l + \eu i u l )-f = 
c cr c 4 



D[t[U[i] ,U[1]] [x]+l/c"2 p[U[i]] [x]u[U[l]] [x] +l/c~2 

u[U[i]] [x]q[U[l]] [x]+l/c"4 e[x]u[U[i]] [x]u[U[l]] [x] ,x ] - f[U[i]] 

Balance of spin 

8i (scs ijl + ^spa ij u l + \u [i spb> ]l + ^spv^u 1 ) - (t^ + ^p [i u 3] + ^u [i q j] 
\ c z c 1 c 4 J \ cr c z 



D[scs[U[i] ,U[j] ,U[1]] [x] + 

l/c~2 spa[U[i] ,U[j]] u[U[l]] + 1/(2 c~2)u[U[i]] spb[U[j] ,U[1]] [x] - 
1/(2 c~2) u[U[j]] spb[U[i] ,U[1]] [x] + 
1/(2 c~4) u[U[i]] spv[U[j]] u[U[l]] - 
1/(2 c~4) u[U[j]] spv[U[i]] u[U[l]],x] - 
1/2 Ct[U[i],U[j]] + l/c~2 p[U[i]] u[U[j]] + l/c~2 u [U [i] ] q[U[j] ] - 
t[U[j] ,U[i]] + l/c~2 p[U[j]] u[U[i]] + l/c~2 u[U[j]]q[U[i]] ) 

Balance of entropy 

di(su l + -se l ) > 



DCs [x]u[U[l]] [x] + l/c'2 se[U[l]] [x] ,x] 
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2.3 Balances on State Space 

As the constitutive functions are defined on the state space, the derivatives have to be performed 
by use of the chain rule. The equations are "labeled" for further use. 

8 l (nu l ) = 



TB = D[nO[x] u[U[l]] [x] ,x] 



ImpB=\ [Kappa] 

D[t[U[i] ,U[1]] [zO[x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] + 
l/c~2 p[U[i]] [x]u[U[l]] [x] + 

l/c~2 u[U[i]] [x]q[U[l]] [zO [x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] + 
l/c-4 e[x]u[U[i]][x]u[U[l]][x],x] - f[U[i]] 

di ( scs ijl (Z) + ^-u [i spb> ]l (Z) + -4« [< spuV J - ( t [iJ] + -kp [i u j] + ^u [l q j] 
\ c 2 c 4 / \ c 2 c 2 



SpinB=D[scs[U[i] ,U[j] ,U[1]] [zO [x] ,zl [x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] + 
l/c-2 spa[U[i] ,U[j]] u[U[l]] + 

1/(2 c~2)u[U[i]] spb[U[j] ,U[1]] [zO[x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] - 
1/(2 c"2) u[U[j]] spb[U[i] ,U[1]] [zO[x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] + 
1/(2 c~4) u[U[i]] spv[U[j]] u[U[l]] - 
1/(2 c~4) u[U[j]] spv[U[i]] u[U[l]],x] - 
1/2 (t[U[i] ,U[j]] + 1/C2 p[U[i]] u[U[j]] + l/c~2 u [U [i] ] q[U[j] ] - 
t[U[j],U[i]] + 1/C2 p[U[j]] u[U[i]] + l/c~2 u[U[j]]q[U[i]]) 



EntropieB= 

D[s [zO[x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5 [x] ]u[U[l] ] [x] + 
l/c~2 se[U[l]] [zO[x] ,zl[x] ,z2[x] ,z3[x] ,z4[x] ,z5[x]] ,x] 
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2.4 Matrix Formulation 



Now the components of the matrices for the matrix formulation of the balances have to be 
collected. 

Matrix of higher derivatives 



xo 


=z0 ' [x] 


XI 


=zl ' [x] 


X2 


=z2 ' [x] 


X3 


=z3 ' [x] 


X4 


=z4 ' [x] 


X5 


=z5 ' [x] 



Then the higher derivatives have to be set correctly by hand. 

(zO) ' [x_] :=nO[L[l]] 
(zl)' [x_] :=u[U[n]] [L[l]] 
(z2)' [x_] :=p[U[n]] [L [1] ] 
(z3)' [x_] :=e[L[l]] 
(z4) ' [x_] :=spa[U[n] ,U[m]] [L[l]] 
(z5) ' [x_] :=spv[U[n]] [L [1] ] 
(nO)' [x_] :=nO[L[l]] 
(u[U[i]])'[x_] :=u[U[i]] [L[l]] 
(u[U[j]])'[x_] :=u[U[j]] [L[l]] 
(u[U[l]]) J [x_] :=u[U[l]] [L[l]] 
(p[U[i]]) J [x_] :=p[U[i]] [L[l]] 
(p[U[j]])' [x_] :=p[U[j]] [L[l]] 
(e)'[x_] :=e[L[l]] 

Then the coefficient matrix 

A00=Coeff icient [TB,zO J [x]] 

AOKJoefficientDrB^l' [x]]+g[L[n] ,U[1]] Coefficient [TB , (u [U [1] ] ) ' [x]] 

A02=Coeff icient [TB,z2 J [x]] 

A03=Coeff icient [TB,z3 J [x]] 

A04=Coeff icient [TB,z4 J [x]] 

A05=Coeff icient [TB,z5 J [x]] 

A10=Coef f icient [ImpB^O' [x]] 

All=Coefficient[ImpB,zl , [x]]+g[L[n] ,U[i]] Coeff icient [ImpB, (u[U[i]] ) ' [x]] + 

g[L[n] ,U[1]] Coefficient [ImpB, (u[U[l]]) ' [x]] 
A12=Coeff icient [^6^2' [x]]+g[L[n] ,U[i]] Coeff icient [ImpB, (p[U[i]] ) ' [x]] 
A13=Coef f icient [ImpB, z3' [x]] 
A14=Coeff icient [ImpB,z4' [x]] 
A15=Coef f icient [ImpB, z5' [x]] 
A20=Coeff icient [SpinB,zO J [x]] 

A21=Coeff icient [SpinB,zl J [x]]+g[L[n] ,U[i]] Coefficient [SpinB, (u[U[i]] ) ' [x]] + 

g[L[n] ,U[j]] Coefficient [SpinB, (u[U[j]]) ' [x]] + 

g[L[n] ,U[1]] Coefficient [SpinB, (u[U[l]])' [x] ] 
A22=Coeff icient [SpinB, z2' [x]]+g[L[n] ,U[i]] Coefficient [SpinB, (p[U[i]] ) ' [x]] + 
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g[L[n] ,U[j]] Coefficient [SpinB, (p[U[j]] ) ' [x] ] 
A23=Coeff icient [SpinB ,z3' [x]] 
A24=Coeff icient [SpinB, z4' [x]] 
A25=Coeff icient [SpinB, z5' [x]] 

Coefficient matrix of dessipation inequality 

\ [Alpha] 0=Coeff icient [EntropieB, zO ' [x] ] 
\ [Alpha] 1= 

Coefficient [EntropieB, zl ' [x]]+ 

g[L[n] ,U[1]] Coefficient [EntropieB, (u[U[l]]) ' [x]] 
\ [Alpha] 2=Coeff icient [EntropieB, z2' [x]] 
\ [Alpha] 3=Coeff icient [EntropieB ,z3 ' [x] ] 
\ [Alpha] 4=Coeff icient [EntropieB , z4 ' [x] ] 
\ [Alpha] 5=Coeff icient [EntropieB ,z5' [x] ] 

The residual matrix 

BO=Simplify[ 

TB - (AOOzO'M + A01 zl ; [x]+A02 z2'[x]+A03 z3'[x]+A04 z4'[x] + 
A05 z5'[x])] 
Bl=Simplify[ 

ImpB - (AlOzO'W + All zl'[x]+A12 z2'[x]+A13 z3'[x]+A14 z4' [x] + 
A15 z5'[x])] 
B2= Simplify [ 

SpinB - (A20z0 , [x] + A21 zl'[x]+A22 z2'[x]+A23 z3'[x]+A24 z4'[x] + 
A25 z5'[x])] 



\[Beta]= Simplify [ 
EntropieB - ( 

\ [Alpha] zO'[x] + \ [Alpha] 1 zl ' [x] +\ [Alpha] 2 z2 J [x] +\ [Alpha] 3 z3'[x] + 
\ [Alpha] 4 z4' [x] +\ [Alpha] 5 z5'[x])] 
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2.5 Liu Relations 



Now Liu's proposition (based on Farkash's lemma) can be applied. The result are the Liu 
equations 

liul= \ [Lambda] A00 + \ [Lambda] 1 A10 + \ [Lambda] 2 A20 -\ [Alpha] 

liu2= \ [Lambda] A01 + \ [Lambda] 1 All + \ [Lambda] 2 A21 -\ [Alpha] 1 

liu3= \ [Lambda] A02 + \ [Lambda] 1 A12 + \ [Lambda] 2 A22 -\ [Alpha] 2 

liu4= \ [Lambda] A03 + \ [Lambda] 1 A13 + \ [Lambda] 2 A23 -\ [Alpha] 3 

liu5= \ [Lambda] A04 + \ [Lambda] 1 A14 + \ [Lambda] 2 A24 -\ [Alpha] 4 

liu6= \ [Lambda] A05 + \ [Lambda] 1 A15 + \ [Lambda] 2 A25 -\ [Alpha] 5 

and the residual inequality 

\ [Lambda] BO + \ [Lambda] 1 Bl + \ [Lambda] 2 B2 -\[Beta] 

The Lagrange parameters can be determined 

lsg=Solve [{liul==0 , liu2==0 , liu5==0} , {\ [Lambda] , \ [Lambda] 1 , \ [Lambda] 2}] ; 
10b=\ [Lambda] /.lsg[[l]] 
llb=\ [Lambda] 1 /.lsg[[l]] 
12b=\ [Lambda] 2 /.lsg[[l]] 
llb=\ [Lambda] 1/ . lsg [ [1] ] ; 
\ [Lambda] 0=10b 
\ [Lambda] 1=1 lb 
\ [Lambda] 2=12b 

Now one can insert the Lagrange parameters into the Liu equations and one gets the restric- 
tions to the constitutive equations. The residual inequality determines the entropy production 
density. 
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